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FOREWORD 

The  mathematical  techniques  presented  in  this  paper  were 
developed  for  the  Defense  Atomic  Support  Agency  as  part  of  the  work 
performed  under  interagency  work  order  804-69.     This  paper  deals 
with  the  specific  task  of  estimating  arbitrary  altitude  variations  in  the 
ion  content  of  the  atmosphere  in  the  nuclear  environment  and  of  develop- 
ing mathematical  techniques  for  calculating  ground  waves  in  such  an 
environment. 

The  listing  of  the  FORTRAN  computer  program  is  available 
from  the  author  upon  request. 
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ABSTRACT 

The  atmosphere  can  apparently  be  ionized  to  such  a  degree 
that  the  ground  wave  is  disturbed,   as  a  result  of,   for  example,   a 
nuclear  event,  and  ionizing  radiation  of  sufficient  intensity  to  overcome 
the  ion -ion  recombination  processes  may  cause  the  ions  to  accumulate 
near  the  surface  of  the  ground.     Gamma  radiation  from  the  fallout  of  a 
nuclear  weapon  may  produce  such  ionization  for  a  considerable  time 
after  the  nuclear  event.     The  problem  is  similar  to  that  of  terrestrial 
radio  wave  propagation  via  the  ionosphere,   except  that  in  the  former 
case  the  ionization  has  a  maximum  near  the  surface  of  the  ground  and 
decreases  exponentially  with  altitude.     In  other  words,   we  have  an  up- 
side down  ionosphere  to  consider.    A  new  theoretical  approach  to  such 
ground  wave  problems  based  on  methods  of  numerical  analysis  has 
been  developed  for  computer  application.     The  computer  computation 
concept  is  discussed. and  sample  computations  based  on  the  fallout 
model  profile  are  presented.    A  particular  fallout  model  based  on  data 
given  in  the  U.  S.  Atomic  Energy  Commission"  s  handbook  on  "The 
Effects  of  Nuclear  Weapons"  is  considered  in  some  detail,   and  the 
effect  of  an  ionization  profile  (ion  number  density  as  a  function  of 
altitude)  on  the  propagation  of  the  ground  wave  is  ascertained.    A 
physically  possible  but  not  likely  extra  attenuation  of  12  dB  at  100  km 
distance  is  noted. 
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GROUND  WAVE  PROPAGATION  IN  AN  IONIZED  ATMOSPHERE  WITH  ARBITRARY 
VARIATION  OF  THE  CONDUCTIVITY  WITH  ALTITUDE 

J.  Ralph  Johler 

1.     INTRODUCTION 

In  an  earlier  paper  (Johler,    1969a)  ground  wave  propagation  in 
a  uniformly  ionized  atmosphere  was  considered  in  some  detail.     It  was 
found  that  uniform  ionization  of  the  atmosphere  in  the  vertical  direction 
could  affect  ground  wave  propagation  if  the  ionization  were  sufficiently 
intense.     Some  nuclear  events  apparently  provide  sources  which  could 
cause  a  disturbance  of  the  ground  wave.     In  this  paper,   we  take  a  more 
realistic  and  more  complicated  approach  by  introducing  a  profile  or 
variation  of  the  atmospheric  conductivity  with  altitude.     For  example, 
gamma  radiation  from  the  fallout  of  a  nuclear  weapon  may  produce 
intense  ionization  near  the  surface  of  the  ground  for  a  considerable 
time  after  the  nuclear  event.     The  ionization  will  decrease  rapidly  with 
altitude.     This  problem  is  similar  to  that  of  terrestrial  radio  wave 
propagation  via  the  ionosphere,   but  in  the  latter  case  the  ionization  has 
a  maximum  near  the  surface  of  the  ground  and  decreases  approximately 
exponentially  with  altitude.     Thus  fallout  produces  an  upside-down 
ionosphere. 

The  basic  unit  for  measuring  fallout  (AEC,    1962,   p. 3 7 5)  is  the 
roentgen>which  is  the  quantity  or  "dose"  of  gamma  radiation  or  x-rays 
that  will  give  rise  to  2.  08(109)  ion  pairs/  cc  under  standard  atmospheric 
conditions  at  sea  level.    An  interesting  case  is  cited  by  AEC  (1962, 
p. 462),   which  gives  the  radiation  level  for  the  BRAVO  test  explosion 
(15  Mt  yield).     At  a  distance  of  190  km  and  a  time  of  96  hours  after  the 
explosion,    the  level  of  3000  roentgens  was  noted.     The  roentgen  can  be 
converted  into  an  average  production  rate, 


roentgens    _    nnt*n9,         roentgens    „    non.9,  .         ,       ,  .,    ,. 

q  =  £^| 2.08(109)  =     98  -  |6Q0    2.08(l09)ions/cc/sec.  (1.1) 


The  unit  dose  rate  may  be  greater  than  the  average  does  rate.     Refer- 
ring again  to  AEC  (fig.    920,   p. 423),   we  can  conclude  that  the  production 
rate    q    should  be  increased  by  a  factor  of  96/  8.  2,   or  approximately  12. 
Thus  the  unit  time  dose  rate  deduced  from  3000  roentgens  at  96  hours 
is  3000/  8.  2,    or  366  roentgens  per  hour   at  one  1  hour  and  1.  5  roentgens 
per  hour  at  96  hours.     The  average  dose  rate  is,   of  course,   between 
these  values  or  3000/  96  or  31.  25  roentgens  per  hour.     (For  references 
of  this  type  see  AEC,   fig.    9.  16a).     The  profile  in  figure  1, which  is 
lbased  on  figure  9.131  in  the  AEC  handbook,   gives  the  dose  rate  variation 
with  altitude  of  a  uniform  distribution  of  radiation  sources  over  an 
infinite  plane.     Thus  a  profile  of  ion  density  and  conductivity  as  a 
function  of  altitude  can  be  deduced  from  the  formulas 


**■-<[$ 


N.  =       -?■      ,  (1.2) 


and 

a  =  2-^ N±_    i2±_t£^     .  (1.3) 

m+     V+j0  e0  me 


Here,     e2  /  m+     is  the  ratio  of  electronic  charge  squared  to  the  mass  of 
the  positive  ion.     The  permittivity  is  given  by    e0  ~  8.  85419(10"     ),  and 
the  positive  ion-neutral  collision  frequency    v+j0     is  given  in  figure  1- 
The  collision  frequency    V+i0  is  multiplied  by  the  mass  ratio  of  neutral 
mass  m      to  the  sum  of  positive  ion  and  neutral  mass  to  account 
for  the  redistribution  of  velocity  after  collision  (Johler  and    Berry,   1965). 
This  point  is  often  overlooked  because,   in  the   case  of  electron-neutral 
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Theoretical  profile  of  atmospheric  conductivity, 
positive  ion  density* and  collision  frequency.     Ion 
densities  are  based  on  fallout  measurements 
270  km  from  ground  zero  and  96  hours  after  the 
BRAVO  nuclear  test.     This  is  the  case    Fn  =  1, 
with  an  average  dose  rate  and  a  nominal  ion -ion 
recombination  coefficient      a 
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collisions,   the  mass  ratio  times  the  collision  frequency  in  the  hydro  dynamic 

equation  is 


m,  +  m0      V*'  °    ~    V*'  °     ' 


since 


Thus, 


and 


m#  <  <  m0    . 


m0  <—  5.  313(10"  6)    for  oxygen, 


m.  -  9.  108(1 0"31)   in  kilograms. 
On  the  other  hand, 


m+  +  m0     +  »•        2      +» °    ' 

if  we  assume  positive  ions  and  neutrals  have  the  same  mass.     Thus, 
strictly  speaking,   the  effective  collision  frequency  is  half  the  actual 
collision  frequency  in  the  case  of  ion-ion  collisions  in  which  each  of  the 
colliding  particles  have  the  same  mass.     This  means  that  the  drift 
velocity  change  after  a  collision  is  redistributed  equally  between  the 
two  colliding  particles.     Since  the  ion  collision  frequency  is  now  more 
precisely  known,   the  mass  ratio  described  above  should  be  taken  into 

account  in  the  equation  of  motion  of  the  ions,  as  was  done  by  Johler  and 
Berry  (1965). 

The  conductivity    cj    reflects  the  effect  of  the  mass  ratio  on  the 
resultant  ion  velocity  change.     The  number  density    N+     is  critically 

3    / 

dependent  upon  the  value  of  the  ion -ion  recombination  rate,     at   cm    /  sec, 
near  the  ground.     Considerable  uncertainty  as  to  the  value  of    at     seems 
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to  exist  in  the  literature.     The  lowest  possible  value  is  the  two -body- 
ion -ion  recombination  rate  used  at  high  altitude,     ot.j  ~  3(1 0~8)  (Knapp, 
19^9,    private  communication,   G.E.    Tempo,   Santa  Barbara,    California); 
the  highest  value  quoted  for  the  three-body  ipn-ion  recombina- 
tion rate  is    a±  ~  3(10"   ).     This  difference,   by  a  factor  of  100,  between 
these  values  means  that  the  profile  in  figure  1  would  change  from  one 
that  would  have  little  effect  on  the  propagation  of  the  ground  wave  to 
one  that  would  seriously  alter  it.     Until  recently  the  three -body  ion-ion 
recombination  rate  most  quoted  in  the  literature    at  ~  2.  65(10"   )  (see, 
for  example,    Loeb,    1961)  was  apparently  based  on  researches  by 
Thomson  (1906,    1924),   Thomson  and  Thomson  (1933),   Rutherford  (1913) 
and  Natanson  (1959).     Experimental  verification  of  the  theory  by  Sayers 
(1938)  gave  the  value  quoted  above,   which  agreed  with  theory „     More 
recently,   Bruechner  (1964)  quotes  a  value  of    3„  95(10"7  )cm3  /  sec  for 
at  .     This  is  based  on  a  corrected  theory  of  ion-ion  recombination  that 
accounts  for  an  error  in  the  Thorns  on -Natanson  theory  by  retaining  the 
more  rigorous  expression  for  the  value  of  the  average  energy  in  the 
analysis.     Brueckner  concludes:     "The  disagreement  (with  the  Thomson- 
Natanson  theory)  may  be  due  to  the  simplifying  assumptions  of  the 
Thomson  theory  as  modified  in  this  paper  to  include  dissociative  col- 
lisions,   such  as  the  failure  to  consider  the  internal  degrees  of  freedom 
of  the  ions  and  the  dispersion  in  energy  and  momenta  of  the  colliding 
particles.     The  changes  which  result  from  a  more  complete  theory  than 
previously  used  suggest,   however,    that  the  previous  agreement  between 
theory  and  experiment  has  been  fortuitous".     More  recently  experi- 
mental evidence  has  been  obtained  in  support  of  Brueckner'  s  theory 
(Hirsh,    1967),   indicating  that  a  considerably  lower  value  of    ai     should 
be  used  near  the  ground. 

Taking  into  account  the  uncertainty  in    04     and  the  production 
rate    q^we  can  translate  the  profile  shown  in  figure  1   simply  by  using 
an  ion  number  density  multiplication  factor    Fn.       Thus  a  change  in  the 
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production  rate  or  the  ion-ion  recombination  rate  merely  modifies  the 
number  density    N+     in  (1.2).     Physically  possible  values,  according 
to  this  author'  s  estimates,   range  from    Fn  =  0.  1    to    Fn  =  20    in  a 
nuclear  environment.     This  would  correspond  to  ground  level  production 
rates    q/q0   =  0.  01  to  400,     where    q0     is  the  average  production  rate 
for  the  BRAVO  profile  (fig.    1).     Thus,   in  the  lower  limit,     Fn  =  0.  1  , 
the  techniques  described  in  this  paper  recovered  the  numerical  values 
calculated  from  the  classical  ground  wave  theory  within  a  fraction  of  a 
decibel  in  the  amplitude  at  a  frequency  of  100  kHz  and  distances  of 
100  km  or  less.     The  particular  case  of  BRAVO  given  by  AEC  would  be 
represented  approximately  by    Fn   =  1  or  2,    depending  on  the  time. 
Therefore,    (1.2)  becomes 


N+   =Fn    /  ■£*.      ,  (1.4) 

where 

qe   =  1.8(107)    f  (1.5) 

and 

afcQ   =  3(10"8)     .  (1.6) 

This  paper  solves  the  electromagnetic  theory  problem  of  a 
ground  wave  signal  propagating  from  an  electrical  point  source  dipole 
to  distances  less  than  approximately  300  km  and  over  finitely  conduct- 
ing ground  with  an  arbitrary  variation  with  altitude  of  the  conductivity 
of  the  atmosphere.     Numerical  values  for    Fn    between  0.  1  and  40  are 
considered  at  a  frequency  of  100  kHz. 
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2.     THEORETICAL  CONSIDERATIONS 

The  problem  of  ground  wave  propagation  at  medium  and  low 
frequencies  in  an  ionized  atmosphere  ordinarily  concerns  short  distances, 
say,    50  to  100  km.     It  is  therefore  more  practical  from  a  theoretical 
point  of  view  to  regard  the  earth  as  flat  and  solve  the  problem  in  terms 
of  a  Cartesian  coordinate  system. 

An  ionized  atmosphere  can  be  accommodated  in  the  analysis  of 
the  ground  wave  with  a  system  of  boundary  conditions  as  depicted  in 
figure  2.     Thus,   the  profile  in  figure  1  is  approximated  as  a  function  of 
altitude  with  a  set  of  plasma  slabs.     This  is  the  method  used  by  Johler 
and  Harper  (1962).     At  the  low  altitudes  that  pertain  to  the  ground  wave, 
atmospheric  reflections  are  a  consequence  of  the  presence  of  ions  that 
exhibit  isotropic  behavior  at  medium  and  low  frequencies.     Therefore, 
advantage  can  be  taken  of  the  simplification  resulting  from  zero 
terrestrial  magnetic  field.     A  detailed  treatment  of  the  effects  of  ions 
is  discussed  in  an  earlier  paper  (Johler  and  Berry  1965). 

The  scheme  for  calculating  the  reflection  process  is  depicted  in 
figure  2.     Each  slab  of  the  system  of  slabs  contains  a  reflection  coeffi- 
cient   Ti     at  the  top,   a  reflection  coefficient    Rj     at  the  bottom,   and  a 
transmission  coefficient    Ui     connecting  it  to  the  adjacent  slab.    As  the 
thickness  of  each  slab    Zi     is  decreased  and  the  number    p    of  slabs  is 
increased,    convergence  can  be  obtained  on  a  reflection  coefficient    T# ,  . 
The  phase  of  the  reflection  coefficient    T.#     is  referenced  to  the  height 
h.     This  reflection  coefficient  is  a  function  of  the  slab  reflection  and 
transmission  coefficients    Ti  ,  Rt  ,  Ut  ,     and    T    within  the  system  of 
slabs.     The  first  slab  of  thickness    h    contains  the  source  dipole  and 
hence  must  be  treated  separately  to  accommodate  the  spherical  waves 
emanating  from  such  an  electrical  point  source.     The  analysis  pre- 
sented below  then  relates  the  bottom  of  the  slab  system  at  height    h 
(fig.    2)  to  the  ground  and  to  a  source  dipole  (transmitter). 
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Figure  2.      Theoretical  treatment  of  the  profile  based  on  the  concept 
of  continuous  stratification  of  the  altitude  variation  of  a 
atmospheric  ion  content  or  conductivity. 
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The  reflection  coefficient    T,,     implies  vertically  polarized 
waves  characteristic  of  a  vertical  electric  dipole  source.     A  similar 
reflection  coefficient    TBm,     for  horizontal  dipole  sources  could  also  be 
used.     Finally,    since  an  isotropic  ionosphere  is  assumed,    the  abnormal 
components  given  by  Johler  and  Harper  (1962)  are  zero  (T.B  =  TB  ,  =  0). 

If    8  (fig.    2)    or    <pj     is  an  angle  of  incidence  for  a  single  semi- 
infinite  slab,   then 


rHP^ 


cos  e  -  ^-l  I  i 

T"   " .  (2.1) 


cose  +  ^i.(^sine)2 

Similarily,   the  ground  reflection  coefficient  for  vertical  polarization 
can  be  written 


k-i     /  ,      f  k.! 


,so     .  '--!      /    i       i    pi-    sinO 
R     = k-2    V  Vk-2     J—     .  (2.2) 


•♦feVl-(fe*0 


The  wave  number  of  the  ground,   assumed  to  be  of  infinite  extent  in  the 
negative  z -direction  (fig.    2)  can  be  written 


k-2   =-Je2    -  i  -1-       ,  (2.  3) 


c   V  e0uo 


where     e2     and    °    are  £ke  dielectric  constant  and  conductivity  of  the 
ground,    respectively.     The  frequency    f    is  given  by    f  =  uo/2TT  (hertz). 
The  speed  of  light  is    c  -  2.  997925(l08)m/  s,   and    e0     is  the  permittivity 
of  space,     e0   =  8.  (10  12).     In  a  particular  slab  along  the  profile, 
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,      _  u>     /  ,       .    o 
1   "  c   V   X   "  1  7>      •  (2.4) 

where 

a  N»  e2  ,  N+  e3  N.  es 

1       *n,(v;>0  +  iui)       m+K^+iuu)      nx.(v:)0  +  iuj)      '     U*  ^ 

/       _         m0  _ 

^&   o  -  _i_  »     a  -  e,  -r  ,    -     , 

rru   +  m„ 


and  where  the  ratio  of  charge  squared  of  the  electron,    e,   to  m  (m, ,  or 
for  the  ion    m^  _.) ,   is  given  by 


e2 
-^-=  2.  81  78(1 0-8)         (for  electrons)     , 


and 


e2 
=  4.  83 10(10  _12)        (for  ionized  oxygen) 


Thus,   the  number  densities    N,  ,  N+  ,     and    N_     depend  upon  the  ionizing 
sources.     At  the  very  low  altitudes  that  pertain  to  the  ground  wave,   the 
electrons  can  be  neglected,   and  the  positive  ions  (profile  illustrated  in 
fig.    1)  provide  the  important  contribution  to  the  finite  conductivity  of 
the  atmosphere.     The  appropriate  collision  frequency  as  a  function  of 
altitude  for  use  in  the  conductivity  calculation  is  also  given  in  figure  1 
as  a  function  of  altitude. 

Since  the  bottom  slab  of  wave  number    k_x     contains  the  source 
of  spherical  waves,   the  propagation  in  this  slab  is  not  as  simple, 
Sommerfeld  (1909,    1926)  and  Weyl  (1919),  and  more  recently  Brekhovskikh 
(i960),   quickly  recognized  that  the  ground  wave  emanating  as  a  spherical 
wave  from  an  electrical  point  source  can  be  expressed  as  a  sum  of  plane 
homogeneous  waves  and  plane  inhomogeneous  waves  that  could  be  con- 
veniently used  to  satisfy  the  boundary  condition  for  the  spherical  wave 
both  at  the  source  and  at  the  plane  boundaries  illustrated  in  figure  2. 
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Thus  it  is  not  necessary  here  to  make  approximations  for  the  boundary 
conditions  at  the  source  as  is  often  done  in  mode  theory  (see  Budden, 
1961,  for  an  extensive  treatment  of  excitation  factors  for  modes).    Based 
on  Sommerfeld1  s  (1909)  basic  approach  all  waves  emanating  from  the 
dipole  can  be  included  in  the  analysis.     This  is  especially  important  in 
dealing  with  the  problem  of  the  ionized  atmosphere  or  where  the 
ionosphere  is  depressed  abnormally. 

For  convenience,   let    (3  =  k_x  ,     the  wave  number  of  the  lowest 
plasma  slab  that  contains  the  source  and  the  observer,   as  shown  in 
figure  2.     It  is  well  known  that  the  vertical  electric  field    Ez  ,     for 
example,   can  be  calculated  from  the  Hertz  vector  for  a  vertical  dipole 
source  given  by 


« °;  =[  iwl°1^^  exp(;ipD)  z  (2.6) 

z       L      -4tt  i  (3      J  D 


by  the  operation 

E*  =[p2+!?Jn"  •  (2-7) 


where 


D2   =  x2  +  y2  +  z2   =  r2  +  z2    ,  (2.  8) 


and  where    I0  I    is  the  dipole  current  moment  in  ampere -meters,   and 
l-i0     is  the  permeability,     \lQ    -  4tt(10~7)  henry/ m,   and    z    is  the  unit 
vector  in  the  z -direction. 

The  factor    exp(-i|3D)/D    can  be  written  as  a  Fourier  integral 
to  show  the  equivalent  plane  wave  structure  of  the  spherical  wave 
source: 
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sEi^£L=^n^-i(Bxx+B7y+B,B,]3i^J  (2.9) 

-00      -00 

where 

3Z  =[p2  -  02  -3y2]"    >  (2.10) 

and 

3X  =  p  sin  0  cos  cp      , 

(3y  =  P  sin  9  sincp       , 

3Z  =  p  cos  6    , 

where     9    and    <p    are  the  vertical  and  horizontal  spherical  coordinate 
angles  respectively.     Bessel'  s  function  of  argument    z    and  order    0  can 
be  written 


2TT 

J0(z)=y      I     exp[±iz  costcp-cpj)]  dcp  •  (2.11) 

o 

Neglecting  the  reflections  from  regions  above  the  source,    (T##  =  0    in 
fig.    2),   one  then  finds  after  considerable  ado  the  integral  form  of  the 
ground  wave  (see,   for  example,   Brekhovskikh  i960,   pp.    242-244  for 
integrals  of  this  type)  to  be 

Ez   =[-f  I«^k^f]     I         Jo(3rsin0){exp[-ip(z-zo)  cos9] 


+  ^e(0)exp[-i0(z+zo)cos9]  j  sin39  d9       .     (z  £  z0)  (2.12) 

The  integration  contour  is  taken  from    0    to    tt/  2    and  thence  from    tt/  2 
to  the  infinity  of  the  imaginary  axis.     It  is  apparent  that  the  solution  can 
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be  split  into  two  integrals:     (a)  the  finite  integral  from    0    to    n/  2    comp^ 
rising  the  homogeneous  plane  waves  and    (b)  the  infinite  integral  from 
n/  2    to    io°     comprising  the  inhomogeneous  plane  waves.     Thus, 


TT 


Ez   =[-32Io^^]     J     {exp[-ip(Z-z0)cos9] 


+  7?e  exp  [-ip(z  +  z0)  cos  9]}  Jo(0r  sin  9)  sin3  9  d9 

00 

+  [-ip2Io^l^j    J  (exp[-3(Z-z0)  sinh9] 

o 
+  R[  exp  [-0(z+zo)sinh9]  \  J0($rcosh9)   cosh39d9  (z  £  z0) 

(2.13) 

where,   in  the  second  integral, 

-isuxhe-l^i-^coshe)* 

7J,'(9)  = ===-         •  (2-14) 

■J1-(bCMh!)' 


-i  sinh9  +  r^ 


The  height  of  the  transmitter    z0     or  the  receiver    z    above  ground 
appears  in  the  exponential  functions,   where  it  is  assumed    z  ^  z0  .     The 
second  integral  is  the  most  difficult  to  evaluate,but  it  converges  prim- 
arily as  a  consequence  of  the  exponential  functions.     In  fact,   the  greater 
the  altitude  of  the  receiver,   the  more  rapid  the  convergence.     As  the 
receiver  altitude  is  reduced  to  zero,   if  the  transmitter  altitude  is 
assumed  to  be  zero,   the  computation  depends  upon  the  Bessel  function 
J0 (z)    for  convergence. 
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If  we  assume  that  the  complete  reflection  coefficient    T,,     at  the 
top  of  the  bottom  slab  in  figure  1  has  been  evaluated,     T,  e  =  F(T,  Tx  ,Ri 
Ui ,  T2  ,R2  ,U2    •  •  •  ).    .Assigning  the  appropriate  ion  density  to    (3  =  k_x 
(thus    (3    may  be  a  complex  number),   we  can  obtain  the  complete  solution 
to  the  problem  by  analogy  to  ground  wave  theory  (2.  13): 

n 
2 

Ez  =\_-flo*<^-~\    I    {  exP  L-iP(z-Zo)  cos  6] 

o 

+  E.  exp[-i$(z  +  z0)cos9]  +  T,,  exp[ -ip(2h  -  z  -  z0)  cos  9] 

+  Bm  T,.   exp[-i3(2h-z+z0)cos6]]-  |l  -7?.  T,  .  exp(-2i  (3h  cos  9)} 

X  J0    Or  sin  9)  sin3  9  d9 

00 

+  [-ip2I0-t^]    \    |exp[-3(z-z0)sinh9] 

o 

+  R[  exp[-0(z+  z0)  sinh  9]  +  T,'.  exp  [  -£(2h-z-z0)  sinh  9] 

i 
+  El  Tj,  exp[-p(2h-z+  z0)  sinh9]}{l   -  J?.'  Ti.exp(-2ph  sinh 9)] 

XJ0(P    r  cosh  9)    cosh39  d9  ,  (2.15) 

(z*  z0) 


where    RB    and  i?^    are  given  by  (2.  2)  and  (2.  14)  respectively.     For  a 
single  semi-infinite  slab,     Tee     is  given  by  (2.  l),and  T,  .     can  be 
written 


-l 
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-isinh  0    -  -=i- 
ki       ' 

A-Of 

cosh  8   j 

-i  sinh  8  +  — — 

h-Ct 

\2 
cosh  8  j 

T.\  = .  (2.16) 

-i  sinh  9  +  7=1.  7  1    -  f  ^=i-  cosh  8  }" 
ki      W  V    ki  y 

The  second  integral  representing  the  inhomogeneous  waves 
converges  for  the  same  reasons  as  in  the  case  of  the  ground  wave  (2.  13). 
In  the  complex  z -plane,     J0  (z)     grows  without  bounds  and  may  cause  the 
second  integral  (2.15)  to  diverge  if    Im  z    is  sufficiently  great  to  over- 
come the  exponential  functions.     This  requires  a  shift  of  the  integration 
contour  to  keep  the  argument    z    of    J0  (z)    real.     Thus,   if 

z  cosh  |  =  (x  +  iy)  cosh  8  =  A  exp(icp)  cosh  8  ,  (2.  17) 

the  required  integration  contour  of  the  second  integral  (2.  15)  in  its 
complex    5  =  5     +  i  ?       plane  is  given  by 

%'  +  iS*  =  cosh"1   [exp  (-icp)  cosh  8}     ,  (2.18) 

5'  +  i  I"  =  In  [exp  (-icp)  cosh. 8  +  J  exp  (-2  icp)  cosh  2  6  -1  }  ,      (2.19) 

provided  that 

■,/c'    ,    .-//»       r        exp  (-icp)  sinh  8  ~l    ._ 

d(§     +!§)=[    -"      F        T/  .,  j  d8     .  (2.20) 

J  exp(-2  icp)  cosh2  8  -1 

The  integration  contour  then  passes  from  zero  to  the  point    P,    in 
figure  3.     At    P,  , 


>i  +i6x'  , 
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Complex 
0-Plane 
0  =  0'+i0" 


Figure  3.      Map  of  complex    8  -plane,   illustrating  integration  contour. 


where  we  may  choose 


B"i  +  iBi    =  In  [exp(-icp)  +  Vexp(-2icp)  -  1  ] 


-1 


=    cosh      [exp(-icp)]  . 
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The  contour  path  is  then  from    Pi     to    P3  ,     where,     at    P3  , 


63  +  i93   =  £n[exp(-icp)  cosh  9  +  J  exp(-2  icp)  cosh2  9  -1],  and    8     is  a 

TT 

select  value    ^  —  .     The  contour  must  now  be  connected  from    P3     to    P2  , 
and    P2     in  the  complex    6     +  i  9       plane  is  given  by 

02  +  iej  ~-\  +  !(§'  +  i5tf)  . 

or 

and  fi 

e£  =  £'   , 

provided    9     in  (2.  19)  is  equal  to  zero. 

Finally,    the  contour  extends  to  the  infinity  of  the  imaginary  axis  along 

,         TT 
a  contour     9    <  —     .     The  break  in  the  contour  between    Pi     and    P2     can 
2 

be  taken  at  a  convenient  point    P3     corresponding  to  a  zero  of  the  Bessel 
function    J0  (|3  r  sin  90),    so  that    90   ^  9   .     Note  that  the  region  between 

TT 

the  original  contour    0  to  —   +  i<»    and  this  new  contour  is  analytic. 
Alternatively,   the  points    Px     and    P3     can  be  projected  to  the     9    axis, 
and  the  contour  from    0    to    P0  ,     P0     to     P2  ,     and    P2     to     9     =  i»    can 
also  be  used.     Finally,  a  contour  is  followed  a  short  distance    0    to    P4 
along  the  real  axis,   then  along  a  straight  line  to    Pi  ,     and  along  an 
identical  contour  to    P3  ,     P2  ,     and  the  infinity  of  the  imaginary  axis. 
All  integration  contours  tested  with  the  computer  program  gave  identical 
results  within  specified  accuracy. 

It  has  been  demonstrated  (Johler,    1969c)  that  integrals  of  the 
type  (2.  15)  involving  Bessel  functions  of  real  argument  (or  with  relatively 
small  imaginary  parts  in  the  argument)  can  be  evaluated  as  an  alternating 
series  in  which  each  term  is  calculated  from  a  Gaussian  quadrature.     The 
zeroes  of  the  Bessel  function  determine  the  limits  of  integration  and  the 
turning  point  for  the  sign  change  of  the  resultant  alternating  series. 
Since  this  technique  is  the  same  as  the  one  discussed  in  an  earlier  paper 
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(Johler,    1969c)  it  will  not  be  repeated  here.    A  complete  computer 
program  for  (2.  15)  is  discussed  in   the  appendix.    This  program  uses  the 
technique  of  continuous  stratification  illustrated  in  figure  1.     The 
computer  subroutine  for  calculating    T#  .     was,   of  course,   required  to 
accommodate  complex  angles  of  incidence    <pt   =  6.     The  theoretical 
basis  of  this  subroutine  is  discussed  elsewhere  (Johler  1967).     The 
complexity  of  the  terrestrial  magnetic  field  has,   obviously,been  re- 
moved from  the  program,    since  this  analysis  is  concerned  with  an 
isotropic  plasma. 

3.     DISCUSSION  AND  CONCLUSIONS 

The  ground  wave  was  calculated  from  (2.  12)  and  compared  with 
the  classical  ground  wave   computer  calculation  (Johler,   1969b;  Johler  and 
Morganstern,    1965)  over  both  plane  and  spherical  ground.     At  distances 
greater  than  about  200  km  the  spherical  correction  may  of  course  be- 
come significant.     However,   at  the  distances  we  are  concerned  with 
here  the  plane  ground  error  was  a  fraction  of  a  decibel  in  the  amplitude. 
The  computer  program  described  in  section  5  which  evaluates  (2.  13) 
was  then  used  in  conjunction  with  the  profile  given  in  figure  1,   together 
with  the  independent  variable    Fn.     Note  that  the  production  rate  is 
related  to  the  square  of  this  factor,     q  «  F|,     according  to  (1.2).     Values 
of   Fn    from  0.  1  to  40,   which  in  this  author1  s  estimate  represent  all 
physically  possible  values  resulting  from  fallout,   were  then  introduced. 

The  most  notworthy  feature  of  the  results  obtained  in  figure  4 
is  a  decrease  in  the  amplitude  of  the  propagated  wave  relative  to  the 
normal  ground  wave  for  values  of    Ffi     of  interest  to  this  problem.     This 
is  physically  a  consequence  of  the  fact  that  the  propagated  energy  is 
slightly  absorbed  by  the  ions  near  the  surface  of  the  ground.     If  the 
atmosphere  were  uniformly  ionized,  the  wave  would  be  much  more 
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severely  absorbed  as    Fn    was  increased  (Johler  1969b).     Since  the 
atmosphere  contains  only  a  thin  film  of  ionization  near  the  ground,    the 
wave  propagates  through  the  air  but  is  concentrated  near  the  ground. 
This  is  true  at  least  for  values  of    FB    up  to  20.     At  this  point,   the 
attenuation  has  increased  to       12  dB  relative  to  the  normal  ground 
waves,   which  would  be  equivalent  to  increasing  the  production  rate  by  a 
factor  of  400  with  respect  to  the  original  BRAVO  model,     Fn  =  1.     The 
possibility  of  such  extreme  production  rates  does  not  seem  likely,   but 
if  it  were  possible  to  increase  the  number  density  of  the  BRAVO  profile 
to    Fn  =  40  there  would  be  a  tendency  toward  excitation  of  a  trapped 
wave  in  the  ionized  atmosphere.     This  would  cause  an  enhancement  of 
the  field.     Also,   the  attenuation  of  the  normal  ground  wave  as  a  function 
of  distance  at  a  steady  12  dB  /  100  km  would  be  replaced  by  a  standing 
wave  as  a  function  of  distance.     It  would  therefore  seem  that  the  12  dB 
attenuation  at  a  distance  of  100  km  is  a  reasonable  maximum  attenuation 
for  physically  possible  values  of    Fn     applied  to  the  BRAVO  model. 
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APPENDIX 
Computer  Program 

The  computer  program  GIDSPAC3  shown  in  the  computation  plan, 
figure  5,   calculates  the  vertical  electric  field    E    by  means  of  (2.  15). 
The  program  has  two  independent  variables,   the  distance  DIST  and  the 
profile  multiplier    Fn   =  PROMLT.     The  amplitude     JEJ   V/m    in  decibels, 
the  phase    Arg  (E)     in  radians,   and  the  phase  correction  in  radians, 

cpc   = r  -  Arg(E)  (radians)  (5»  1) 

are  given,   where  the  field,     E,     is  defined  by 

E  =  |e|   exp[iuot  -  i  -  r  -icpc]  .  (5.  2) 

Homogeneous  waves  EZHA,  EZHP,   PHICl  and  inhomogeneous  waves 
EZIHA,   EZIHP,   PHIC2  are  given  separately  for  convenience,   together 
with  the  total  field,  EZA,   EZP,   PHIC3.     Appropriate  comments  explain 
each  step  in  the  program.    A  feature  of  the  program  is  the  shift  in  the 
integration  contour  shown  in  figure  3,   which  tests  the  precision  of  the 
numerical  analysis  and  interpolation  routines  used.     The  shift  can  be 
controlled.     Thus,   in     (2.17),  <p  =  AK1P,   and  with  AK1P=    0, 
the  contour    0,  P4,P0,P2     to    i«    is  followed.     With  AKlP  =  AKlP, 
the  contour    0,  P4 ,  Px ,  P3  ,   P2     to    i«    is  followed.     Another  example 
would  be  AKlP  =  AKlP/  2,   which  would  follow  a  contour  between  the 
above  contours  up  to  the  point    P2  .     The  contour  from    P2     to    i®    is 
never  disturbed  since  it  has  been  optimized  according  to  (2.  18). 
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Figure  5.       Conceptual  plan  of  program  GIDSPAC3. 
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